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SUMMARY 

The assumed natural strain (ANS) formulation of finite elements has undergone rapid de- 
velopment over the past five years. The key formulation step is the replacement, in the po- 
tential energy principle, of selected displacement-related strains by independently assumed 
strain fields in element natural coordinates. These strains are not generally derivable from 
displacements. This procedure was conceived as one of several competing methods to solve 
the element locking problem. Its most noteworthy feature is that, unlike many forms of 
reduced integration, it produces no rank deficiency ; furthermore, it is easily extendible to 
geometrically nonlinear problems. Many original formulations were not based on a varia- 
tional principle. The objective of Part I is to study the ANS formulation from a variational 
standpoint. This study is based on two hybrid extensions of the Reissner-type functional 
that uses strains and displacements as independent fields. One of the forms is a genuine 
variational principle that contains an independent boundary traction field, whereas the 
other one represents a restricted variational principle. Two procedures for element-level 
elimination of the strain field are discussed, and one of them shown to be equivalent to the 
inclusion of incompatible displacement modes. In Part II, the 4- node C° plate bending 
quadrilateral element is used to illustrate applications of this theory. 
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1. INTRODUCTION 


The ass um ed natural strain (ANS) formulation of finite elements is a relatively new de- 
velopment. A restricted form of the method was introduced in 1969 by Wiliam [14], who 
constructed a 4-node plane-stress element by assuming a constant shear strain indepen- 
dently of the direct strains and using a strain-displacement mixed variational principle. 
A different approach advocated by Ashwell [l] and coworkers regarded “strain elements” 
as a way to obtain appropriate displacement fields by integration of assumed compatible 
strain fields. These and other forms of assumed-strain techniques were overshadowed in 
the 1970s by developments in reduced and selective integration methods, but have recently 
begun to attract attention [2,6,8,10,13]. The primary motivation behind recent work has 
been the construction of simple and efficient finite elements for plates and shells that are 
locking-free, rank sufficient and distortion insensitive, yield accurate answers for coarse 
meshes, fit naturally into displacement-based programs, and can be easily extended to 
nonlinear and dyn ami c problems. Elements that attain these attributes are collectively 
known as high performance elements. 

Over the past 20 years investigators have resorted to many ingenious devices to construct 
high-performance elements. Among the most successful ones we can mention patch-test- 
verified incompatible displacement models, reduced and selective integration, mixed and 
hybrid formulations, stress projectors, the free formulation, and assumed natural strains. 
The underlying theme is that although the final product may look like a standard displace- 
ment model so as to fit naturally into existing finite element programs, the conventional 
displacement formulation is abandoned. (By “conventional” we mean the use of conforming 
displacement assumptions into the total potential energy principle.) 

Another common historic trend is that certain deviations from the conventional formulation 
were initially made without variational justification and in fact labelled as “variational 
crimes” by applied mathematicians. In some cases such as reduced numerical integration, 
reconciliation was achieved later after surprisingly good results prompted explanation. 
In other cases, notably non-conforming elements and the patch test, a comprehensive 
mathematical theory is still in the making. 

The present paper seeks to interpret the assumed natural strain (ANS) formulation from 
a variational standpoint. The justification is based on hybrid extensions of the Reissner- 
type functional that uses the strains and displacements as independent fields. We restrict 
our considerations to linear elasticity although the straightforward extension to geometric 
nonlinearities is one of the strengths of the ANS formulation. In Part II, the 4-node C 
plate-bending quadrilateral is used as a specific example to illustrate the application of the 
present variational interpretation. 
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2. PROBLEM DESCRIPTION 


2.1 Governing Equations 

Consider a linearly elastic body under static loading that occupies the volume V . The body 
is bounded by the surface S t which is decomposed into S : S u U SV Displacements are 
prescribed on S u whereas surface tractions are prescribed on St. The outward unit normal 
on 5 is denoted by n = n». 

The three unknown volume fields are displacements u = Uj, infinitesimal strains € = <*/» 
and stresses o = o^j . The problem data include: the body force field f = /» in V , prescribed 
displacements u = u » on S ut and prescribed surface tractions t = t» on Sf 

The relations between the vol um e fields are the strain-displacement equations 


c = |(Vu + V T u) = Du or 

= JK/ + «/,<) 

(1) 

(where superscript T denotes transposition), the constitutive equations 


o = E c or Oij 

= Eiikiiki »n V, 

(2) 

and the equilibrium (balance) equations 



—div o = D *o = f or 

o%j\] ■(■ /i = 0 in V, 

(3) 


in which D* = — div (divergence) denotes the adjoint operator of the symmetric gradient 
D = §(V + V r ). 

On S the surface stress vector is defined as 

o n = o.n, or o n i = Oijnj. (4) 

With this definition the traction boundary conditions may be stated as 

o n = t or a ij n 3 = U on St* (5) 

and the displacement boundary conditions as 

u = u or « » = u» on S u . (6) 

2.2 Notational Conventions 

An independently varied field will be identified by a letter without superscript, for example 
u, c, o. A dependent field is identified by writing the independent field symbol as super- 
script. For example, if the displacements are independently varied, the derived strain and 
stress fields are denoted by 

_ i(v + V T )u = Du, ff “=E<“ = EDu. (7) 

Given a finite element subdivision of V , quantities pertaining to the e th element will be 
identified by superscript (e), for example u^ e ^, wherever appropriate. At an interface 
between two elements e and /, superscripts (ef) and (/e) will identify interface quantities 
considered as part of e and /, respectively. 
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3. THE HU-WASHIZU AND REISSNER FUNCTIONALS 


In the conventional Hu-Washizu functional the displacements u, stresses a and strains e 
are independently varied. Arranging the strain and stress components as vectors, and the 
elastic moduli in E as a matrix, the functional may be expressed as| 

L(u,c,<r) = J [{( T E£ + f T (t’‘-()-f T u] dV ~J s dS ~J s * TndS - 

From L one obtains the conventional stress-displacement Hellinger-Reissner functional by 
eliminating € through the inverse of (2), namely € = €° = E 1 o . Another Reissner- 
type, strain-displacement functional is obtained by eliminating o through the constitutive 
relation (2), namely o = v* = E«, which yields 

R( u, £)=/[- i£ T E£ + £ t E<“ -fn]dV- jT (<) J- (u -H)dS- £ t T u iS. (9) 
Setting c = c u reduces R to the potential energy functional 

P(u) = J [i(c“) T E£“ - f r u] SV - Kf (u -A)dS- J s t T u dS, (10) 
generalized with a S u term over its usual expression. 


4. HYBRID FUNCTIONALS 
4.1 Independent Boundary Tractions 

If the functional (9) is used to construct finite elements, the displacement field u should 
be continuous in V because of the presence of £ u , whereas the assumed strain field may 
be discontinuous. To account rigorously for displacement discontinuities it is necessary to 
add the interelement surface tractions t as new independent field which plays the role of 
Lagrange multiplier. Let Si denote the union of interelement boundaries traversed twice 
(one for each adjacent element); on Si neither displacements nor tractions are prescribed. 
Then R expands to the hybrid functional 

tf(u,£,t) = R(u,£) - / t T u dS. (11) 

JSi 


f There are several equivalent statements of this functional, differing from one another in 
transformations based on the divergence theorem. For example in Gurtin [5, p. 122] the 
stress divergence appears. Some authors attribute this specific functional to B. Fraeijs de 
Veubeke, who indeed published a version of it in 1951, four years before Hu and Washizu. 
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Figure 1. Sample finite element mesh to illustrate 
computation of integrals in H 


For later reference we note the specialization € s of (11) to the generalized potential 
energy functional of Jones [12] 


P(u,t) = P(u) — [ t T u <fS, 

JSi 


( 12 ) 


where P(u) is given by (10). 


The meaning of the integrals in H may be illustrated on the two-dimensional mesh of 
Figure 1: 

/ =/ +[ + / 

Jv ^ Jv^) Jvw Jvw 7v(») 


X. ? A** + Xi i > + Ssi s) 

I Si = ^ Is}" ~ Jsl 1 ' + /,« + X<*> 

f = £/ =/ + [ +/ + / 

J5 < , J S (ti) ./$(*•»> J $(*•*> 

c,/ • 


(13) 


where element identification conventions stated in Section 2.2 have been followed. It is 
seen that in the integrals over V, 5 U and S t each element appears once, whereas in Si 
adjacent elements appear twice. 
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4-2 First Variation 
The first variation of Hi 

6H = 6„H + 6 t H + S t H , (14) 

♦ 

yields the Euler equations and interelement linking conditions, which are underlined in the 
expressions below. The three components of 6 H are 


S U H = J ( V<x« -jf 6udV + J s « ~ t) r *u dS + j s « - t) T 6u dS , (15) 

6,H = f E(f“ - tf Sc'JV- ( (n-0) r t(Ec)„Ag, (16) 

Jy Js % 

6 t H = [ }£_6tdS. (17) 

JSi 

Note that there are two contributions to the element interface integrals, one from 6 U H 
and another from 6 t H. Putting the parts together and decomposing into element-pair 
contributions we get 


f [« - t) r *u + u T St] dS = Yl f k (e)T ^ (e) ~ 

^ ej Js{tf) (18) 

_t( e '> T $u( e ) - dS . 


In the absence of applied internal tractions, interelement equilibrium requires t (e/) = 
— t^ e ^, which substituted into (16) reduces the right-hand side to 


y; [ \o* M T 6 xiM - c t W T 6 uM - tM T 6{uW - u (/) ) + (u (e) - u* / )) T a‘* / >l dS. 

Ti Js{ " n l ~ 2 

(19) 

If we assume a compatible displacement field, the above equation reduces to 

y / « (e) - < (/, ) T «uW dS t (20) 

T7 ;5(, ’ /) 

which means that the interelement equilibrium condition appears as the Euler equation 
corresponding to the variation of the interface displacements. 

4.3 A Restricted Variational Principle 

If the displacement field b incompatible we should in principle retain t as an independent 
boundary-traction field satisfying t (e/) = over interelement boundaries. One way 

to achieve this is to assume a continuous stress field a* over element boundaries, so that 


tW) = <r*.n (c) = (e) , t (/e) = <r*.n (/) = <r*.(-n (e) ) = -<r* 


(«) 


(21) 
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The presence of an independent boundary traction field is computationally disadvantageous 
because additional degrees of freedom must be retained on elements sides. This contradicts 
one of the tenets of high-performance element construction noted in the Introduction. It 
would be more convenient if o* could be identified with the strain-derived stress field , that 
c * = <r e ss Ec on 5*, because we would have only two independent fields, u and £, as in 
(9). The strain freedoms can be eliminated at the element level as explained in Section 6, 
and we are left with standard displacement connectors. The corresponding functional is 

H(U,€)=R - \ «) T U<*5. (22) 

JSi 

But in general a * n is not continuous between elements. One can argue, however, that con- 
tinuity is achieved in the limit of a converged solution. A variational statement such as 
SH = 0 is called a restricted variational principle [3, Ch. 11] because the governing field 
equations of §2.1 are satisfied only at the exact solution. Away from it, 6 H — 0 gener- 
ally violates interelement-equilibrium field equations although it may provide satisfactory 
numerical approximations. 

Stress-displacement (rather than strain-displacement) functionals of this form have been 
used by Pian and coworkers [11,12], who transform the interface integral into an element 
volume integral and in doing so introduce a stress divergence term. 

4.4 Finite Element Classification 

Finite element models derivable from R , H and H may be classified into several types 
according to the number of independent fields and the continuity conditions on those 
fields. Following are some general comments on the most interesting combinations, which 
are summarized in Table 1. 

1. Continuous displacements. The independent boundary field t is not needed, and 
we can work with the mixed functional R. If the strain field is discontinuous, strain 
freedoms may be eliminated at the element level as explained in Section 6. Continuous 
strains are in principle possible but impractical in general structural applications 
where material interfaces, plasticity, and sudden thickness or area changes may occur. 

2. Discontinuous displacements. The displacement field contains conforming and non- 
conforming portions. Assumed strains axe discontinuous and may be eliminated at the 
element level. Displacement degrees of freedom associated with non-conforming modes 
may be also eliminated if separable. The governing functionals are H or H . With the 
latter an independent traction field t is required; degrees of freedom associated with 
t must be retained at the assembly level. 

In practice elements are often constructed as a combination of these types with conven- 
tional displacement models. Thus part of the strain field may be considered as completely 
derivable from displacements and part as independently assumed, as discussed in Section 
8. This was in fact the scheme originally used by Wiliam [14]. The C° plate bending 
quadrilaterals studied in Part II provide another important example. 
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Table 1. Assumed-Strain Finite Element Models Derivable From R, H and H 


Element 

Type 

Governing 

functional 

Independent 

fields 

Interelement 
continuity on * 
net 

Element 

connected 

fields 

Element 

condensable 

fields 


R 

U,€ 

c 

d 

u 

€ 

(n) 

R 

U,€ 

c 

c 

U,€ 


(m) 

H 

U,€ 

n 


ut 

€ 

(IV) 

H 

U,€,t 

d 

d c 

u^,t 

£ 

* c=continnoufl, d=discontmuou8. 

f conforming part only if separable 

as per (33) 


5. DISCRETIZATION 

5.1 Assumptions 

In this section the finite element discretization of the hybrid functionals H and H is 
studied. That is, we focus attention on element types labelled (HI) and (IV) in Table 1. 
In the sequel it will be assumed that the displacement boundary conditions are identically 
satisfied by u, whence the strain-displacement hybrid functionals reduce to 

H{ u,£,t) = J [c r E(*“ - \e) - f T u] dV - j^fudS- t r u dS. (23) 

ff(u,e) = J ^« t E(«“ — §t) - f r uj dV — J t T udS — J {ff^) T udS. (24) 

The framework used here accomodates both continuous and discontinuous displacements. 
The FE assumption may be written 

u = Nv in V, £ = Aa in V, t = Ts on Si. (25) 

Here matrices N, A and T collect displacement shape functions, assumed natural strain 
functions and interface traction functions, respectively, whereas column vectors v, a and s 
collect nodal displacements, strain amplitudes, and interface tractions amplitudes, respec- 
tively. The derived fields in V are 

6 U = DNv = Bv, <r“ = EBv, = E« = EAa. (26) 
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5.2 Discrete Equations 

On inserting the assumptions (23-24) into (21-22) we obtain the bilinear algebraic forms 

E{v, a,s) = -§a T Ca + a T Pv-v T Ls-v T p, (27) 

H[y, a) = -±a T Ca + a T (P - R)v - v T p = - Ja T Ca + a T Pv - v T p. (28) 

where 


C = f A T EAdV = C T , P= f A t EB dV, L = / N T TdS, 

Jv Jv JSi (29) 

R = f (EA)^N dS, P = P - R, P = J y NTf & + J NT * dS * 

* ft 

Observe that (28) results on substituting Ls by R T a in (27). Making these forms station- 
ary yields the linear systems 


-C 

P T 

0 


p 0 1 f a 1 f 

0 -L ^ V ► = < p ► , 

-L T 0 J ( s j 0 , 

(30) 

P o]{v} = {p}- 

(31) 


for (27) and (28), respectively. In both cases the first matrix equation is the discrete 
analog of (16), and expresses internal compatibility. The second matrix equation is the 
analog of (15) and expresses internal and boundary equilibrium, and, in the case of (31), 
approximate boundary compatibility. The third matrix equation in (30) is the analog of 
(17) and expresses boundary compatibility. 

5.S Displacement Field Decomposition 

With view to further developments the assumed displacement field is decomposed as 


u = u c + u a- (22) 

where u c is continuous (compatible, conforming) in V and uj discontinuous (incompatible, 
non-conforming) on 5,-. It will be further assumed that this decomposition can be effected 
in terms of the shape functions, i.e., 

u = N c v c + N dV d , ( 33 ) 

where the vj freedoms are defined element-by-element and may in principle be condensed 
out. This assumption holds for elements in which non-conforming shape functions are 
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“injected” over a compatible set. For the H functional, as shown in Section 4.2 the S% 
integral exactly vanishes for the conforming displacements: 



(34) 


On the other hand, for E the corresponding 5* integral also vanishes at the converged 
solution. Taking this into account, equations (30-31) expand to 


1 

0 

o 

Cl 

O 

i 


f \ 

a 


' 0 ' 

pj 0 0 0 

i 

V« 

> = i 

P ‘ , 

Pj 0 0 -L d 


Vd 


Pd 

i — 

o 

0 

1 

t- 1 

R.S 

0 

1 


, ■ . 


. o ; 


(35) 


r-c 





(36) 


in which Pd = Pd — Pa. and where c- and d-subscripted matrices and vectors are given 
by integrals similar to (29) in which N is replaced by N c and Nd, respectively. 

6. STRAIN ELIMINATION 

The strain degrees of freedom may be eliminated at the element level by static condensation 
or by enforcing kinematic constraints. These two techniques are studied below. 

6.1 Static Condensation 

This is a well known variationally consistent procedure which will be illustrated for the 
system (30). From the first matrix equation get a at the element level: 

a = C _1 Pv = Q.v. (37) 

Substitution into the second equation gives 

[£• ?]{;}-{:}• m 

where K = P T C -l P = P T Q 4 = Q^CQ, is a stiffness matrix. Similarly, (31) condenses 
to 

Kv = p, (39) 

where K = p T C -1 P = Q^CQ, and Q J = C -1 P. The separable non-conforming degrees 
of freedom Vj, if present, may be condensed out following a similar procedure. 


10 


6.2 Kinematic Constraints 

A second elimination procedure has been used recently in the construction of ANS C° 
plate and shell elements. It will be described by considering the system (35) that displays 
separable conforming and non-conforming displacement shape functions. A kinematic 
constraint that links strain to displacement degrees of freedom is established: 

a = Q„ v c + Qd*d- (40) 

This relation may be constructed by collocation, least-square fitting or some other means. 
Often = 0. For example, in the Bathe-Dvorkin element [2] studied in Part II collocation 
of natural shear strains is done at the quadrilateral midpoints. 

If the following conditions hold: 

(a) the dimension of Va and a are the same so that Pa is square; 

(b) matrix Pa — CQ d is nonsingular; 

then the relation (40) may be interpreted as a variationally-consistent constraint on non- 
conforming displacements. In effect, the first equation of (35) becomes 


(Pc - CQ c ) v c + (Pa - CQa)va = 0, 

(41) 

wn ptitp 

VJ = -(Pj - CQj)-'(P. - CQ> 0 = Wv„ 
a = (Q c + QaW)v c — Qv c . 

(42) 

If (as often happens) Qj = 0, Q = Q c . Replacing the constraints (42) 
form 5’(a,v c , va,t) and setting its first variation to zero yieldsf 

into the discrete 

Ur V]{V}-{".-}- 

(43) 

where 

K* = Q t CQ, L* = W r L d , p* = p„ + W T p d . 

(44) 

Similarly, for (34) we get the stiffness equations 


■£?* 

K v c = p , 

(45) 


where K = Q CQ, in which Q d results on replacing Pa by Pa in (41-42). 

Note that condition (a) above may be relaxed if the dimension of vj exceeds that of a by 
selecting a subset of Vj that satisfies (b), and statically condensing out the remainder. 

f One obtains K* = Q r (2P 0 +2P<*W-CQ) which simplifies to (44) because PjW = CQ— P c . 
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6.8 Relation to the Strain Projection Approach 

If the dimension of a exceeds that of Vj (in particular, if the assumed displacement field 
is conforming) the constraint (40) is in general inconsistent with a strain-displacement 
variational principle. In such a case a connection with other techniques for improving ele- 
ment performance can sometimes be established. For example, suppose that the assumed 
strains £ are constant and equal to ? over each element, and that the displacements are 
continuous. We can choose a = 2, and A = I so that (40) may be written 

l = Bv. (46) 

This is the strain-projection approach, also called averaged-B or the B approach. If B 
is determined by collocation at the element center, (46) is equivalent to one-point re- 
duced/selective integration on the potential energy functional, see e.g. Hughes’s textbook 
[7, Ch. 4]. 

7. LIMITATION PRINCIPLE 

The famous limitation principle of Fraeijs de Veubeke [4] was originally stated for stress- 
displacement mixed finite elements, but holds for many strain-displacement elements as 
well. The principle is applicable when the displacement-derived strain field c u is contained 
in the assumed strain field c: 

€ 3 £ u = Du = Bv. (47) 

This inclusion cam be expressed in matrix form as 

e = Aa = Ba„ + A z a x = [B A,]{*”}. (48) 

Here a„ contains the same number of entries as v whereas A*, which may be empty, 
contains “excess” strain modes. Consider elements of type (III) based on the functional 
JET. Inserting (48) into (30) we get 

-C T „ -C„ c* 

Cw Cj® 0 

O 0 — L r 

where 

C vv = J B t EB dV, C vx = jf B T EA a dV, C xx = jf A^EA dV. (50) 

The first two matrix equations give a„ = v and a x = 0. Hence the system is equivalent to 
(38) in which K = C„„ as simply the potential energy stiffness matrix. Consequently the 
stiffness equations may be directly constructed from the generalized potential energy func- 
tional (12) and the independent strain assumption has no effect. Of course the conclusion 
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only applies if the strain degrees of freedom are solved for in a manner consistent with 
the variational equations (49); for example by static condensation. - If the derived field 
varies over V , assuming a constant strain field Z for £ is a safe way to guard against the 
limitation principle. 

A similar analysis of type (IV) elements on the -ET-derived system (31) shows that the 
limitation principle does not generally hold unlessJEtv = 0 . For arbitrary v this implies 
that the interface integral vanishes, in which case H reduces to the mixed functional R. 


8 . PARTIAL STRAIN ASSUMPTIONS 

It is common practice to assume only part of the strains as independent fields. For example, 
in the C° plate bending element studied in Part II independent assumptions are only 
made for the transverse shear strains whereas the bending strains are entirely derived 
from displacements. The partial strain assumption may be expressed as 

(5i) 

where independent strain assumptions are made only for € a = Aa. For one has £{> = £ 5 . 
The R and H functionals require obvious modification in the volume term; for example, 


R(u,£ a ) 





dV + surface terms 


(52) 


while for H an additional adjustment in the S* integral is required. The resulting principles 
take a particularly simple form if the constitutive coupling term E tt {> and E^ vanish, in 
which case 

S = R B (u,£.) + JMu) ( 53 ) 

where R a is a mixed strain-displacement principle involving £ a , and Pb is a potential-energy 
principle involving the strain energy. 


9. CONCLUSIONS 

The key results of the present study may be summarized as follows. 

1 . The mixed strain-displacement functional of Reissner type, R, can be expanded to two 
hybrid functionals, H and H , to account for incompatible displacements. Whereas 
SR = 0 and 6H = 0 are genuine variational principles, 6H = 0 represents a restricted 
variational principle. 

2. Several types of assumed-strain finite elements may be constructed using R , H or 
H. The most practical elements for inclusion into existing displacement codes are 
those in which (l) strain and non-conforming-displacement degrees of freedom can be 
eliminated at the element level and ( 2 ) avoid surface traction connectors. 
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3. Strain degrees of freedom may be eliminated by static condensation or through kine- 
matic constraints. The latter technique can be presented in a variationally consistent 
form if the conditions stated in Section 6.2 hold, in which case it can be interpreted 
as a constraint on non-conforming displacements. Special versions of this technique 
are closely related to the strain-projection approach. 

4 . DeVeubeke’s limitation principle applies to finite element models derivable from func- 
tionals R and H if the strain elimination procedure is variationally consistent. 

5. The present variational formulations may be readily modified to account for partial 
assumptions on the strain field. 
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A VARIATIONAL JUSTIFICATION OF THE ASSUMED 
NATURAL STRAIN FORMULATION OF FINITE ELEMENTS. 
II. THE C° FOUR NODE PLATE ELEMENT 


CARMELO MILITELLO 
CARLOS A. FELEPPA 

Department of Aerospace Engineering Sciences 
and Center for Space Structures and Controls 
University of Colorado 
Boulder, Colorado 80309-0429, USA 


SUMMARY 

In Part II we use the four-node C° plate bending element to explore some of the possibilities 
opened by the theory presented in Part I. This element is chosen because the version 
presented by Bathe and Dvorkin [2], MITC4 , can be considered the simplest assumed 
natural strain element tha; allows several possibilities to be studied in a straightforward 
manner. We focus our attention on the governing functionals R and H presented in Part 
I, assuming independent s r rain fields only for the transverse shear strains. Besides MITC4 
we consider three formulations (two mixed and one hybrid) that collectively represent a 
variational justification for the assumed strain technique. In addition, we examine reduced 
and selective-integration elements to compare their behavior with that of the present strain- 
assumed elements. 
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1. INTRODUCTION 


1.1 J-Nodc C° Bending Plate Element Formulation 

We start with the formulation of the four-node Reissner-Mindlin plate element whose 
degrees of freedom (d.o.f) are the transverse displacement w and the two rotations 6 X and 
6 V about the x and y axes, respectively, as shown in Figure 1. We expand the displace m ent 
field in the usual way: 

w = N% (r, s) Wi 

6 x =N i {r,s)0 xi (I) 

= Ni{r,s) Qyi 

where 

Ni{r,s) = i(l + r ir )(l+s { s), * = 1,2, 3, 4 (2) 


are bilinear shape functions. The strain field derived from the displacement field is 


vv 

iff* 

l*i 


** 

II 

If.* 

= — z 

$x,y 

-i* 


= w t9 

H 

1 

H 

9 

II 

+ Q v 


( 3 ) 


We take advantage of the decoupling between bending and shear energies by using different 
assumptions for each one. We assume that the bending strains coincide with the bending 
strains computed from the displacement field: 

C XX — °xx 

£ vv = e yv (4) 

= e* v 

The shear strains components in the Cartesian basis x , y, z derived from the displacement 
field are 

= “’•* + IS) 

i;. 

After some manipulations we can obtain the covariant components of the shear strains in 
terms of the natural coordinates r and s as 


Iri = «\r + 0r 

l7* ~ + 0* 


( 6 ) 

( 7 ) 


17 



Figure 1. Element coordinate system and notational conventions 


where 


fir = y,r + *,r 

fi» = -** y,. + *y 


( 8 ) 

(9) 


1.2 The Assumed Covariant Shear Strain 

We consider two different assumptions for the covariant shear strains, 


and 


1 Ism = <*3 


(i - *) - (i + 3 ) 

< 7 r , = <*1 2 — + 2 2 

(1 - f ) ^ „ (i±£) 


+ <* 4 ‘ 


Kra = “1 

ns* - <*2 


( 10 ) 

(n) 


( 12 ) 

(13) 


The bilinear assumption (10)— (11) is of the same form as that proposed in [2j. The constant 
strain assumption (12)-(13) is studied to see whether there are connections to the selective 
reduced integration (SRI) technique discussed by Hughes [3]. 


18 



2. MIXED ELEMENT BASED ON THE FUNCTIONAL R{u,c) 

Up to now we are working with a compatible displacement field and a discontinuous strain 
field. Hence we use the functional R(u, f) presented in §3 of Part I. No boundary field is 
necessary and the constants a,- can be obtained at the element level. 

The element displacement field is 


which can be expressed as 
where 


[ w | 

u = < Ox > , 

UJ 


u = N„v, 


a) bending strains: 


b) shear strains: 



‘ Ni 0 

o 

o 

£ 

• 

• , 
• 

o 

N* = 

0 Ni 

o 

• 

• 

• 

o 

Si 

o 


. 0 0 

s? 

o 

o 

• 

• 

ft 

5? 

*r-(«i **i 

6yi ... U>4 0,4 0y4 ) . 

ed from the displacements are 



’ <- \ 



> = B£v e , 



2e“ j 



h xy / 


T* = 

\ 1 1 ‘ } = b : v c . 



1 J 


(14) 

(15) 

(16) 
(17) 


(18) 


(19) 


The independently varied strains are: 

a) bending strains: the same as obtained from the displacement field, t.e., (18). 

b) shear strains: 

i ={£;}= b > (20) 

Replacing (18), (19) and (20) into the functional R and carrying out the integrations at 
the element level we obtain 


J2(v e ,a) = 


where 


-v T Kf Ci 

2 v c ***6 


Kf c 


' c - la T C“a + vfL“a - vjp° , 

(21) 

/ (BJ) T E»B{<B' 
Jv* 

(22) 
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( 23 ) 


c~ = [ (b;) t e,b; <jv, 

Jv 

L ca = [ (B:) t E,B a 4 dV, 

Jv 

p e = f Njf dV+f Njt dS 

Jv* is; 


(24) 

(25) 


Here vector f collects applied distributed forces conjugate to w , 6 X and 0 V . On performing 
the variations we obtain the matrix equation 


Kr 


■Hil-P) 


(L co ) t -C°° 

From the second equation we obtain the shear strain coefficients 

a = (C“)-'(L“) r Vc = Q c v c 

which replaced into (26) gives the statically condensed system 


(26) 


(27) 


(K r + QfC“Q.)v. = p« (28) 

Here Kj* is the bending stiffness matrix, which is also obtainable from the potential energy 
principle, and QjC aa Q c stands for the new shear stiffness matrix; cf. §8 of Part I. 

Equation (27) can also be obtained by minimizing the following shear energy error norm: 

n. = ± / (i-i”) T E.(i-T,“)<iv 

Jv « 

where vector 7 collects the independent shear strains (lO)-(ll) or (12)-(13), and 7 “ collects 
the shear strains evaluated from the displacement field, equation (19)* The minimization 
of this norm using an independent stress field instead of a strain field was proposed by 
Barlow [ 4 ] as a way of deriving stress-assumed hybrid elements. 

We have implemented two elements based in the form ( 21 ) and the assumptions ( 10 )— (ll) 
and (12)-(13), which will be identified as Pi and Pi, respectively, in the sequel. The 
results obtained for the simple shear and bending tests illustrated in Figures 2 and 3 are 
summarized in Tables 1 and 2 . We have compared these results to those obtained using 
SRI and MIT Cl elements. The results indicate that Pi and Pi behave poorly when 
elements are distorted and that Pi is not equivalent to SRI. 

An interesting result is that if we use one point reduced integration to compute L ca , both 
elements PI and P 4 yield the same results obtained using SRI. 

We can obtain another expression for Q c , called Q* in the sequel, from the field proposed 
by Bathe and Dvorkin [2] for the covariant shear strains. This expression relates four 
strain coefficients a to the nodal degrees of freedom v c . The elements of Q c are given 
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Figure 2. Bending test 



0 X : 0y : 0 


in Appendix A. It is important to realize that Q c obtained for element P 4 matches the 
matrix Q* only for rectangular shapes. Consequently, the variational principle based on 
the functional R justifies the assumed natural strain technique for rectangular shapes. 
However, what can we say about distorted shapes? We need Q c = Q* for all possible 
configurations to generalize that justification. 
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Table 1. Bending Test (FEM/ Theory-Figure 2) 


a 

Node 

MITC4 

w 9 f 

SRI 

w 

PI 

w S f 

P4 

w 9 f 

0 . 

S 

1.00 

1.00 

1.00 

1.00 

1.00 

1.00 

1.00 

1.00 


6 

1.00 

1.00 

1.00 

1.00 

1.00 

1.00 

1.00 

1.00 

1. 

6 

1.00 

1.00 

0.90 

1.00 

0.88 

1.00 

0.44 

0.44 


6 

1.00 

1.00 

1.10 

1.00 

1.07 

1.00 

0.47 

0.47 

2. 

5 

1.00 

1.00 

0.80 

1.00 

0.74 

1.00 

0.23 

0.23 


6 

1.00 

1.00 

1.20 

1.00 

1.06 

1.00 

0.28 

0.29 


Table 2. Shear Test (FEM/Theory-Figure 3) 


a 

Node 

MITC4 

SRI 

PI 

P4 



w 

w 

w 

W 

0 . 

5 

1.00 

1.00 

1.00 

1.00 


6 

1.00 

1.00 

1.00 

1.00 

1 . 

5 

1.00 

1.00 

1.40 

1.00 


6 

1.00 

1.00 

0.85 

1.00 

2. 

5 

1.00 

1.00 

3.06 

1.00 


6 

1.00 

1.00 

0.99 

1.00 


3. INCOMPATIBLE DISPLACEMENTS. THE FUNCTIONAL #(u,£,t) 

Following the general procedure outlined in §6.2 of Part I, we add to the transverse dis- 
placement w the four midside incompatible shape functions of an eight node element. In 
thw way the bending behavior is unchanged. We denote by v d the nodal values associated 
with these “injected” incompatible shape functions. The new displacement field can be 
written as 

«-[N. N„l{£j (29) 

where 

f \{l + r)(l - s 2 ) §(1 - r)(l - s 2 ) §(1 + s)(l - r 2 ) *(1 - s)( 1 - r 2 ) 

oooo 

0 0 0 0 


• (30) 


N d = 
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The bending strains do not change, and for the displacement derived shear strains we have 

'M£H v - +b?v '- <3i) 


If we introduce the new strains into the variational principle, we must use the functional 
J7(u,e,t) because the displacement field will be discontinuous. Then, we have to introduce 
a traction field t over the boundary . This traction field is a (line) shear resultant, and 
for simplicity we shall assume that it is constant on each element side. On performing the 
variations, the following expression at the element level is obtained: 


' K£° 0 p°° L®** 

0 o P°° L* 

(pea)T (p^T -Q** 0 

(I/*) T (L*) T 0 0 


where 

f (B*) t E.B*^ 

Jv 

(33) 


P **= / (B d ) r E*B® dV 
Jv • 

(34) 


L rt = [ N IdS 
Js*t 

(35) 


L* = [ Nj dS 
J S*t 

(36) 


p d = [ Njt dS + f Njf dS 
Js • Jv • 

(37) 

Now imposing the relation 

a = Q:v e 

(38) 

we obtain 

Vd = 

(P da )- r (c aa Q* - (P ca ) T )v c = W c v c . 

(39) 



Replacing both relations in the variational principle and taking variations with respect to 
v c and t, the following expression at the element level is obtained: 


K£ c + Q: t C°°Q: + W C L*] / v c \ _ f p c + Wjp d 1 

(L°* + W c L dt ) T 0 J \ t J \ 0 / 


The stiffness matrix proposed in [2] for the plate element, namely, K£ c + Q* T C aa Q*, 
can be clearly identified in the preceding expression. It is not necessary to compute the 
contribution L** because it comes from the compatible displacement and will cancel with 
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Figure 4. Cantilever beam discretization 


the contribution of the neighboring element. On the other hand, the contribution L from 
the incompatible mode does not vanish. If t vanishes the stiffness matrix reduces to that 
of [2] but the nodal force vector will generally be different. Thus it is worth emphasizing 
that the variational principle gives a consistent treatment for the distributed loads. 

The matrix P^ 1 is singular for rectangular elements, but we know that in this case Q c is 
equal to Q* and there is no need to introduce the incompatible displacement field. 


4. NUMERICAL EXAMPLE 

To check the behavior of the functional IT(u,«,t) we analyze a cantilever beam with 
two distorted elements, as depicted in Figure 4. The assumed independent shear strain 
corresponds to equations (10) and (11). We are interested in two load cases: a uniform 
bending moment at the tip (Figure 2); and a uniform transverse load at the tip (Figure 
~3jr In both cases Poisson’s ratio is set to zero to compare the results to those obtained 
through the Euler-Bemoulli beam theory. 

Uniform Bending Moment . The theoretical solution for this problem requires a linear 
variation for By and a quadratic variation for the transverse displacement xv. The results 
obtained with MITC4 coincide with the theoretical results. So do those obtained with the 
present formulation labeled ANSH (for Assumed Natural Shear Hybrid). 

The value obtained for t is of roundoff error order (1.10“ 12 ). Then, in this case, both 
formulations are equivalent and the work of the incompatibility can be disregarded. 
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Table 3. Normalized displacements (FEM/Theory) for bending, t — 1.25 — 12 


Node 

MITC4 

ANSH 


w 

9, 

W 

9, 

5 

1.000 

1.000 

1.000 

1.000 

6 

1.000 

1.000 

1.000 

1.000 


Table 4. Normalized displacements (FEM/Theory) for shear, t = -2.227 


Node 

MITC4 

ANSH 


w 

9 , 

w 

9s 

5 

0.930 

1.077 

0.892 

1.003 

6 

0.912 

0.920 

0.891 

1.002 


The external load vector is the same for both formulations because the external bending 
moment does not interact with the transverse displacement. 

Uniform Transverse Load. The theoretical solution requires a quadratic variation in 6 V and 
a cubic one in w. In this case we must expect the computed solution to be approximate. 
The results obtained are shown in Table 4. Clearly the ANSH formulation is less sensitive 
to element distortion. The lack of symmetry can be observed at the third decimal position. 
The convergence and symmetry for the rotation is excellent. The value obtained for t is 
not negligible. Note that in this case the external load vector is not the same for the 
MITC4 and ANSH formulations. 


5. CONCLUSIONS 

We have illustrated the theory presented in Part I (lj through the study of several 4-node 

C° plate elements with independently assumed shear strains. The following conclusions 

emerge from this study. 

1. Elements Pi and P4 based on the mixed functional 22(u,c) are variationally impec- 
cable. Pi behaves well in the bending test and P4 passes the shear patch test. Their 
performance deteriorates markedly, however, if the element geometry departs from 
the rectangular one. 

2. The MITC4 element imposes a shear strain- displacement relation (38) obtained by 
midpoint strain collocation. This kinematic relation is not a priori derivable from a 
mixed variational principle such as 8 R = 0. 

3. A variationally consistent modification of MITC4 , named ANSH, is obtained by in- 
troducing incompatible displacement modes and an independent surface traction t (in 
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this case a shear line force), and using the hybrid functional -£f(u,£,t) for the shear 
energy portion. The results are similar to those of MITC4- Although this element is 
more expensive to form, it does provide a consistent treatment of applied distributed 
loads. 

4. The MITC4 element stiffness matrix is recovered by setting the boundary traction 
field t of ANSH to zero. However, the nodal load vector for distributed applied forces 
will generally be different. 

The techniques illustrated here are obviously applicable to the construction of other types 
of strain-assumed elements based on the various functionals presented in Part I [l]. In 
particular, the use of the restricted hybrid principle H, in which the boundary tractions 
are not retained as independent degrees of freedom, remain unexplored. 

A key result of this investigation is that any change in the strain-displacement interpo- 
lation from the variationally consistent interpolation must be associated in some way to 
the addition of incompatible displacement modes. This property is closely linked to the 
limitation principle stated in §7 of Part I. 
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Appendix A 


Bathe and Dvorkin [2] proposed the same kind of shear strain interpolation we have used in equa- 
tions (10)— (11). To determine the coefficients a,- they imposed the following midpoint-collocation 


relations: 



nX + tft 


ai = 

2 ’ 

aa = 


7:. a +7:. s 


a» = 

2 * 

a 4 = 




2 

_+ 

2 


i:, + 


where superscripts 1,2,3, 4 indicate the node where expressions (6) and (7) must be evaluated; 
see Figure 1. Through the application of the relations of Section 1 and after some algebra we 
obtain 


• = q: v« 


where 



Oa Os a* ) 


vf = (tu 0,i 0 9 i ... 0, 4 ) 


• 0 .5 x±--Z2 —0.5 


q: = 


0.5 


l 0 


4 

0 


0 

0 

0.5 


0 

0 


0 

0 


—0.5 *3-^4 


0 

0.5 


hi - y« - *3 

4 4 


- H i *2 - Oa _Q 5 ys -yq yx - . g j 

4 4 4 4 


—0.5 ftir-yi ^-7-^4 

4 4 
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